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1. Introduction 

We introduce and analyze a crude model for the random loose packings of 
granular matter. These packings, as well as random close packings, were carefully 
prepared by Scott et al in the 1960's [S,SK], mainly in samples of steel ball bearings. 
Gently pouring samples of 20,000 to 80,000 spheres into a container, the lowest 
possible volume fraction obtainable - the so-called random loose packing density - 
was determined to be 0.608 ± 0.006. 

The above refers to monodisperse steel spheres immersed in air; similar exper- 
iments were performed with spheres of other materials immersed in other fluids; 
variations in the coefficient of friction and in the effective gravitational force lead 
to somewhat different values for the random loose packing density [SK]. 

Matter is generally described as "granular" if it is composed of a large number 
of noncohesive subunits each of which is sufficiently massive that its gravitational 
energy is much larger than its thermal energy. A common example is sand. 

There are several classic phenomena characteristic of static granular matter, 
in particular dilatancy, random close packing, and random loose packing, none of 
which can yet be considered well-understood; see [dG] for a good review. A basic 
question about these phenomena is whether they are sharply defined or inherently 
vague. Dilatancy has recently been associated with a phase transition measured 
by the response of the material to shear [SN], which answers the question for this 
phenomenon. The case of random close packing is controversial and awaits further 
experiment; see [TT,R]. Our main goal here is to analyze this question with respect 
to random loose packing, to determine whether or not traditional theoretical ap- 
proaches to granular matter predict a sharply defined random loose packing density. 
It is clear that any experimental determination of a random loose packing density 
will vary with physical conditions such as coefficient of friction. It is possible, how- 
ever, that there is a precise geometrical quantity which underlies the phenomenon, 
and it is our goal to investigate this question. 

We will eventually specialize to a certain two dimensional model, but we be- 
gin our discussion more ambitiously. Consider a model using a large collection of 
impenetrable, unit mass, unit diameter spheres in a large container, acted on by 
gravity and with infinite coefficient of friction between themselves and with the 
container. We will not try to understand a random loose packing density as an 
absolute minimum, at which the spheres can form a bulk material. Rather, we will 
use a probabilistic framework in the spirit of Edwards' theory of granular matter 
[EG]. More specifically, we expect, but cannot show, that the following is true, and 
use it as our motivation. Consider a probability density on the set of all mechani- 
cally stable packings of the spheres in their container, with the probability density 
of a packing c proportional to exp[—E{c)], where E{c) is the sum of the heights, 
from the fioor of the container, of the centers of the spheres in the packing c. We 
expect that such an ensemble will exhibit a gradient in the volume fraction, with 
volume fraction decreasing with height, and that there is a well-defined random 
loose packing density as one approaches the top of the packing, where the (ana- 



1 



logue of hydrostatic) pressure goes to zero. More specifically, we expect that as one 
takes an infinite volume limit, the probability distribution for the volume fraction 
of the top layer of the packing becomes concentrated at a single nonzero value. We 
emphasize that we are focusing on a bulk property near the top of the configuration, 
not a surface phenomenon. 

To analyze random loose packing in such a probabilistic framework we make 
several simplifications of the above model, first to an ensemble of those packings 
which are limits, as the gravitational constant goes to zero, of mechanically stable 
packings; we effect this by setting E{c) = in the relative density exp[— £'(c)]. With 
this simplification configurations are now, in their entirety, representative of the top 
layer in the original model. Next we consider the two dimensional version of the 
above: congruent frictional unit disks in mechanically stable configurations under 
vanishingly small gravity. Note that each such disk must be in contact with either 
a pair of supporting disks below it or part of the container. (Here and elsewhere 
in this paper we neglect events of probability zero, such as one sphere perfectly 
balanced on another.) We simplify the model one last time by replacing the disks by 
congruent squares, with edges aligned with the sides of the (rectangular) container, 
each square in contact with either a pair of supporting squares below it or the fioor 
of the container. This is now a granular version of the old model of "(equilibrium) 
hard squares" [H], which is a simplification of "hard disks" and "hard spheres" 
(sec [AH] for a review), in which gravity is neglected but kinetic energy plays a 
significant role. We emphasize that in our granular model there is no longer any 
need to concentrate on the "top layer" ; in fact we will eventually be concerned with 
an infinite volume limit which, as usual, focuses on the middle of the collection of 
squares and lets the boundaries grow to infinity. 

It is this granular model which is the focus of this paper. We have run Markov 
chain Monte Carlo simulations on the model, with the following results. We ini- 
tialize the squares in an allowed configuration of some well-defined volume fraction 
anywhere between 0.5 and 1. If the initial volume fraction is not approximately 
0.76, the packing expands or contracts, relatively quickly, into the range 0.76 ±0.01; 
see Figs. 1 and 2. 

The process is insensitive to the tightness of fit of the initial configuration in 
the containing box except for absolute extremes. If the side walls of the containing 
box abut a closely-packed initial configuration, the simulation cannot significantly 
change the volume fraction; alternatively if the width of the box is much larger 
than that of the initial configuration the simulation will drive the configuration to a 
monolayer on the floor; but both extremes are negligible and the equilibrium volume 
fraction is otherwise insensitive to the fit of the initial configuration in the box. More 
precisely, we found that the equilibrium volume fraction should be accurate if the 
floor length is between 2\/N and SVn, where N > 100 is the number of squares. 
Since we will be conjecturing the behavior of the model in the inflnitc volume limit, 
the equilibrium conflguration should be a single bulk pile, so the floor length should 
be on the order of \/N. To understand the lower bound, note that, at any volume 
fraction, a configuration occupies the least amount of floor space when the squares 
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are arranged in a single full triangle. The bottom level of such a triangle has just 
under \/2N squares. Assume the containing box fits tightly around the triangle. If 
the triangle has volume fraction greater than 0.754 then the configuration will not 
be able to decrease to this equilibrium volume fraction; we avoid this by ensuring 
that the floor length is at least (0.754)~-'^-\/2iV. To arrive at the upper bound we 
performed simulations on fixed particle number and let the floor length vary. We 
saw that the equilibrium volume fraction was reliable so long as the floor length 
was less than about at least for > 100. 

We conclude that, for floor lengths in the aforementioned acceptable range, the 
equilibrium volume fraction is found to depend only on the number of squares in 
the system. The main goal of our work is an analysis of the distribution of volume 
fraction - both the mean and standard deviation - as the number of particles 
increases. We conclude that the limiting standard deviation as particle number 
goes to inflnity is zero, so the model exhibits a sharp value for the random loose 
packing density, which we estimate to be approximately 0.754. 

The heart of our argument is the degree to which we can demonstrate that in 
this model there is a sharp value, approximately 0.754, for the equilibrium volume 
fraction of large systems, and we postpone analysis of error bars to later sections. 
But to understand the value 0.754, consider the following crude estimate of the 
volume in phase space of all allowable packings at flxed volume fraction First 
notice that the conditions of the model prevent the possibility of any "holes" in 
a conflguration. Furthermore, if we consider any rectangle in the interior of a 
conflguration, each horizontal row in the rectangle contains the same number of 
squares. (One consequence is that in the inflnite volume limit each individual 
conflguration must have a sharply deflned volume fraction; of course this is quite 
distinct from the degree of spread of the volume fraction among all conflgurations.) 
Now consider a very symmetrical conflguration of squares at any desired volume 
fraction 0, with the squares in each horizontal row equally spaced, and gaps between 
squares each of size (! — (/))/(/) centered over squares in the next lower horizontal row; 
see Fig. 3. Consider these squares to represent average positions, fix all but one 
square in such a position, and consider the (horizontal) degree of motion allowed 
to the remaining square. There are two constraints on its movement: the gap 
size separating it from its two neighbors in its horizontal row, and the length to 
which its top edge and bottom edge intersects the squares in the horizontal rows 
above and below it. These two constraints are to opposite effect: increasing the 
gap size decreases the necessary support in the rows above and below. A simple 
calculation shows that the square has optimum allowed motion when the gap size is 
1/3, corresponding to a volume fraction of 0.75, roughly as found in the simulations. 
In other words, this crude estimate suggests that the volume in phase space (which 
we effectively estimate, for N squares, to be where L is the allowed degree of 
motion of one square considered above) is maximized among allowed packings of 
fixed volume fraction by the packings of volume fraction about 0.75. 

To obtain accurate physical measurements a method of sedimentation has been 
developed to prepare samples of millions of grains in a controlled manner; see [JS] 
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and references therein for the current state of the experimental data. In these 
experiments monodisperse grains sediment in a fluid. The sediment is of uniform 
volume fraction, at or above 0.55 depending on various experimental parameters. 
To achieve the low value (0.55 ± 0.001 [JS]), the grains need to have a high friction 
coefficient and the fluid needs to have mass density only slightly lower than the 
grains, to minimize the destabilizing effect of gravity. (In the absence of gravity 
one could still produce a granular bed by pressure; we do not know of experiments 
reporting a random loose packing value for such an environment.) 

Our results suggest that whatever the initial local volume fraction of the bed, 
on sedimentation (in low effective gravity) most samples would have a well-defined 
volume fraction, the random loose packing density of about 0.55, with no intrinsic 
lower bound on the sharpness of the value. This should be verifiable by a sequence 
of experiments of increasing dimensions. 

There have been previous probabilistic interpretations of the random loose 
packing density, for instance [MP,PC] . A distinguishing feature of our results is our 
analysis of the degree of sharpness of the basic notion, which, as we shall see below, 
requires unusual care in the treatment of error analysis. 

2. Results 

We performed Markov chain Monte Carlo simulations on our granular model, 
which we now describe more precisely. We begin with a fixed number of unit edge 
squares contained in a large rectangular box B. A collection of squares is "allowed" 
if they do not overlap with positive area, their edges are parallel to those of the box 
B, and the lower edge of each square intersects either the fioor of the box B or the 
upper edge of each of two other squares; see Fig. 3. Note that although the squares 
have continuous translational degrees of freedom in the horizontal direction, this 
is not in evidence in the vertical direction because of the stability condition: the 
squares inevitably appear at discrete horizontal "levels" . 

Markov chain simulations were performed as follows. In the rectangular con- 
tainer B a fixed number of squares are introduced in a simple "crystalline" con- 
figuration: squares are arranged equally spaced in horizontal rows, the spacing 
determined by a preassigned volume fraction </>, and with squares centered above 
the centers of the gaps in the row below it; see Fig. 4. The basic step in the simula- 
tion is the following. A square is chosen at random from the ciirrent configuration 
and all possible positions are determined to which it may be relocated and produce 
an allowed configuration. Note that if the chosen square supports a square above it 
then it can only be allowed a relatively small horizontal motion; otherwise it may 
be placed atop some pair of squares, or the floor. So the boundary of the conflgu- 
ration plays a crucial role in the ability of the chain to change the volume fraction. 
In any case the positions to which the chosen square may be moved constitute a 
union of intervals. A random point is selected from this union of intervals and the 
square is moved. The random movement of a random square is the basic element 
of the Markov chain. It is easy to see that this protocol is transitive and satisfles 
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detailed balance, so the chain has the desired uniform probability distribution as 
its asymptotic state [NB]. See Fig. 5 for a configuration of 399 squares after 10^ 
moves. Our interest is in random loose packing, which occurs in the top (bulk) layer 
of a granular pile, and we assume that the entirety of each of our configurations 
represents this top layer. We emphasize that our protocol is not particularly appro- 
priate for studying other questions such as the statistical shape of the boundary of 
a granular pile, or properties associated with high volume fraction, such as random 
close packing. 

After a prescribed number of moves, a volume fraction is computed for the col- 
lection of squares as follows. Within horizontal level Lj, where j = corresponds 
to the squares resting on the floor, the distances between the centers of neighboring 
squares is computed. (Such a distance is 1 + g where g is the gap between the 
squares.) Suppose that nj of these neighboring distances are each less than 2, and 
that the sum of these distances in the level is Sj. At this point our procedure will be 
complicated by the desire to obtain information during the simulation about inho- 
mogeneities in the collection, for later use in analyzing the approach to equilibrium. 
For this purpose we introduce a new parameter, p. For fixed < p < 1 we consider 
those levels, beginning from j = 0, for which rij is at least 0.75p times the length of 
the box's fioor. Suppose Lj^p^ is the highest level such that it, and all levels below 
it, satisfy the condition. We then assign the volume fraction 

to the assembly of squares. (The factor 0.75 represents the volume fraction we 
expect the box's fioor to reach in equilibrium. Note that any two such calculations 
of volume fraction of the same configuration may only differ by a term proportional 
to the length of the boundary of the configuration, so any inhomogeneity is limited 
to this size.) Such a calculation of volume fraction was performed regularly, after 
approximately 10^ moves, producing a time series of volume fractions (/»t for the 
given number of squares. (We suppress reference to the variable p for ease of 
reading. As will be seen later all our results correspond to the choice p = 0.4, so 
one can, without much loss, ignore other possibile values.) Variables and (j)t+i 
are highly dependent, but we can be guaranteed that if the series is long enough 
then the sample mean: 



1 ^ 

t=i 

will be a good approximation to the true mean of the target (uniform) probability 
distribution for the given number of squares [KV]. 

We created such time series (/)t, each of about 10^ terms (roughly 10^° elemen- 
tary moves), using values p = 0.2, 0.4, 0.6, 0.8, on systems for the following numbers 
of squares: 100m, and 1000m, for m = 1,2, ...,9, with varying initial volume 
fractions. For each system we needed to determine the initialization period - the 
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number of moves necessary to reach equilibrium - and then the total number of 
moves to be performed. Both of these determinations were made based on variants 
of the (sample) autocorrelation function f(k) of the time series {(pt | 1 < ^ < T} of 
volume fractions, defined for 1 < /c < T by: 

I T-k+l 

/w^ T-fe+i ^ {4>t-mt+k-^), (3) 

where is the mean of the series. This function is easily seen to give less reliable 
results as k increases, because of limited data, so it is usual to work with functions 
made from it as follows. One way to avoid difficulties is to restrict the domain of 
/; we define the "unbiased" autocorrelation fi{k) by fi{k) = f{k) for k < T/10. 
Another variant we consider is the "biased" autocorrelation /2, defined for all 1 < 
A; < T by: 

. T-fc+l 

h{k) = j, E {<t>t-mi+k-4>). (4) 

which reduces the value of f{k) for large k. (See pages 321-324 of Priestley [P] for 
a discussion of this biased variant.) We consider both variants of autocorrelation; 
to refer to either we use the term fj. 

With these autocorrelations we determined the smallest k = kz such that 
fj{k) = 0. We then computed the sample standard deviation af^ away from zero of 
fj restricted to k > kz, and defined kj to be the smallest k such that \fj{k) \ < af.; 
see Fig. 6. (For ease of reading we sometimes do not add reference to j to quantities 
derived using fj.) This defined the initialization period. Then starting from (f)j.j 
we recomputed the autocorrelation fj and aj, and determined the mixing time, the 

smallest k = kM such that \fj{k)\ < aj . kM was interpreted as the separation k 
needed such that the random variables and (pt+k are roughly independent for all 
t > k. (We performed the above using the different definitions of volume fraction 
corresponding to different values of p, allowing us to analyze different geometrical 
regions of the samples. For each system of squares we selected, for initialization 
and mixing times, the largest obtained as above corresponding to the various values 
of p, which was always that for p = 0.2, corresponding to the lowest layers of the 
configuration.) 

Once we determined kj and kM we ran the series to (pp, where F = kj + Tku 
for some T > 20. The values of kj and kM are given in Tables 1 and 2; the empirical 
means and standard deviations of volume fraction are given in Tables 3 and 4. 

In all our results we use p = 0.4 to minimize the boundary effects presumably 
associated with small or large p. (With large p the lowest level may have undue 
influence on the volume fraction; with small p, the surface levels could have undue 
influence. Note that the arrangements of squares on the lowest level and the surface 
levels are not restricted by the arrangements of squares below and above them, 
respectively, and so the corresponding volume fractions are not bound to the logic. 
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discussed above, which suggested that each level should equilibrate at a volume 
fraction of about .75. In spite of this, we found that using any series corresponding 
to p in the range 0.2 < p < 0.8 generated a similar result.) For all systems the 
volume fraction quickly settles to the range 0.76 ± 0.01 and we can easily see from 
Table 4 that the empirical standard deviations decrease with increase of particle 
number. In Fig. 7 we plot the empirical standard deviations against particle number, 
and in Fig. 8 the data is replotted using logarithmic scales. In Fig. 8 we include 
the best least-squares fit to a straight line y = ax + obtaining a = —0.5004 and 
b = —0.8052. The corresponding curve is included in Fig. 7. Also included in both 
graphs are 90% confidence intervals for the true standard deviations, obtained as 
described in the next section, using /2. (There was not enough data to obtain a 
confidence interval by this method for the system with 8995 squares.) The same 
data is reanalyzed in Figs. 9 and 10 with confidence intervals derived using /i. 

We use the close fit to the line in Figs. 8 or 10, corresponding to 33 data points 
in a range of particle number varying from 100 to 9000, to extend the agreement to 
arbitrarily large particle number, and therefore to claim that the standard deviation 
is zero in the infinite volume limit, or that there is a sharp value for the random loose 
packing density. The argument is supported from the theoretical side by noting the 
closeness of the slopes in Figs. 8 and 10 to —1/2. A slope of —1/2 would be expected 
if it were true that an equilibrium configuration of N squares could be partitioned 
into similar subblocks which are roughly independent - a proposition which would 
not be surprising given a phase interpretation of granular media [R] . Verifying such 
independence might be of some independent interest but would require much more 
data and much longer running times. 

This is our main result, since it shows how to make sense of a perfectly well- 
defined random loose packing density within a granular model of the standard 
Edwards' form. 

As to the actual asymptotic value of the volume fraction in the limit of large 
systems, we assume that our simulations sufi'er from a surface error proportional to 
-\/]V for a system of N squares. The least-squares fit of a function of form A+B/y/N 
to the data (see Figs. 11 and 12) yields A — 0.7541, and the good fit suggests an 
(asymptotic) random loose packing density in our granular model of about 0.754. 

Our argument concerning asymptotically large systems depends on the fit of 
our standard deviation data to a curve, and the degree to which this fit is convincing 
depends on the confidence intervals associated with our simulations. In the next 
section we explain how we arrived at our confidence intervals. 

3. Simulation data analysis 

A good source for common ways to analyze the data in Markov chain Monte 
Carlo simulation is chapter 3 in Newman and Barkcma [NB]. We will give a more 
detailed analysis, following the paper by Geyer [Gey] in the series put together for 
this purpose by the statistics community [Gel]. As will be seen, our argument is 
based on the precision of estimates of various statistical quantities, and necessitates 
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a delicate treatment. 

Our simulations produce a time series Cj of (dependent) random configurations 
of squares. From this we produce other series g{cj) using functions g on the space 
of possible configurations c, in particular the volume fraction gi{c) — (f) and g2{c) = 
(0 — K)"^ for constant K. 

We use the common method of batch means. As described in the previous 
section, we first determine an initialization time kj and a mixing time /cm for our 
series Cj from autocorrelations. After removing the initialization portion of the 
series, we break up the remaining W terms of the series into w > 2 equal size 
consecutive batches (subintervals) , each of the same length W/w, discarding the 
last few terms from the series if w does not divide W evenly. 

It should be emphasized that rarely, if ever, are conclusions drawn from a finite 
number of Monte Carlo simulations a literal proof of anything interesting. We are 
going to obtain confidence intervals (using the Student's t-test) for the mean and 
standard deviation of the volume fraction of our systems of fixed particle number. 
The t-test 's results would be mathematically rigorous if in our simulations we had 
performed infinitely many moves; of course this is impossible, so we will try to make 
a convincing case that we have enough data to give reliable results. Ultimately, this 
is the most sensitive point in our argument. 

Assume fixed some function gr, and denote the true mean of g{c) by jig. Assume, 
temporarily, that enough moves have been taken for the t-test to be reliable. (We 
will come back to this assumption below.) With the notation g[c) for the empirical 
average (l/iu) Ylik^9ip))k of gic), the variable: 

^-^^ __ (5) 

T.ki{9{c))k-g{.c)Y 

has a t-distribution, allowing one to compute confidence intervals for ^g. 

The above outline explains how (given the validity of the t-test) we could 
compute confidence intervals for the mean value of the volume fraction for the 
time series associated with our simulations for fixed numbers of squares. A small 
variation allows us to give confidence intervals for the standard deviations of these 
variables, as follows. 

Denote the true standard deviation of g{c) by Ug. Using conditioning, 

Prob(/ig e / and Ug & J) = Prob(//j, e /)Prob(crg e J | /i^ G /) 

= Prob(/ig e /) ^Prob((73 e J \ jig e /i)Prob(//g e h \ jig e I), 

i 

where {/j} is a partition of /. We have discussed how to obtain / so that the 
factor Prob(/Ug G /) is at least 0.95. We now want to obtain J so that the factor 
Prob(crg G J I G /) is also at least 0.95, and therefore Prob(//p G / and ag G J) 
is at least (0.95) (0.95) > 0.90. 
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Consider, for each constant K, the random variable 



^K = J{l/w)J2Mc))j-K]^. (7) 



Using (5) with {g{c) — K)'^ playing the role of g{c), we can obtain a 95% confidence 
interval for the mean of E^, which we translate into a 95% confidence interval Jk 
for the mean of E/^. Assume the partition so fine that within the desired precision 
Jk — Ji only depends on i, where K & li. Note that ii K — Hg, then the random 
variable has as its mean the standard deviation ag. So if we let J = J^, then 
Prob((Tp e J \ iJ,g E li) > 0.95 for all z, and therefore Prob(crg e J \ /ig E I) > 0.95. 
In practice the union J = UiJi is easy to compute. 

In the above arguments we have assumed that enough moves have been taken to 
justify the t-test, which has independence and normality assumptions which are not 
strictly satisfied in our situation. We now consider how to deal with this situation. 
Some guidance concerning independence can be obtained from the following toy 
model. 

Assume that for the time series of the simulation one can determine some 
number /cm, perhaps but not necessarily derived as above from the autocorrelation 
f{k), such that variables 0^ and (f^i+k in the time series are roughly independent 
if A; > kM- We model this transition between independent random variables as 
follows. 

Let T and M be nonnegative (integer) constants. For < t < T and 1 < m < 
M — 1 we first define independent, identically distributed random variables X^m 
and from these define: 

XtM+m = (l ~ Jd)^*^ J^-^(t+l)M, (8) 

together defining for < t < TM — 1. Note that variables Xf and Xf+m are 
independent for m > 2M — 1. 

A simple calculation shows that: 

^tM+m = ( ^ ) + ^ )^(t+l)M- (9) 

m=0 

Then another simple calculation shows that: 



1 ^-^"^ 1 ^"^ 1 1 

^^^TM ^ Xm^j;[j2XtM + ^{Xo+XTM) + ^{Xo-XTM)]. (10) 
m=0 t=l 

In other words St is the mean of roughly K independent variables. 

Returning to the question of the assumptions in the t-test, the toy model 
suggests that the independence assumption is easily satisfied. The normality as- 
sumption is usually taken as the more serious [Gey]. But we note from [DL] that 
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the t-test is quite robust with respect to the normahty assumption. Although the 
robustness of the t-test is well known and is generally relied on, in practice one 
still has to pick specific batch partitions in a reliable way. This is not covered in 
[Gey]. We arrived at a standard for batches of length 10 times mixing time for our 
series as follows. In outline, we use mixing times as computed above to standardize 
comparison between our systems with different particle number. Those for which 
our runs constituted at least 800 mixing times are assumed to give accurate values 
for the mean volume fraction. Various initial segments of these runs are then used, 
with various choices of batch partitions, to see which choices (if any) give reliable 
results for confidence intervals. Batches of size 10 mixing times proved reliable even 
for initial segments in the range of 20-100 mixing times, so this choice was then 
used for all systems. We emphasize that we are using this method to determine a 
minimum reliable batch size on the sequence of configurations, and then we apply 
this to the time series (jit as well the time series [cjit — KY. We now give more details. 

For most of the systems of particle numbers 100-900 we have over 500 mixing 
times worth of data, yet for some of the systems of particle numbers 1000-9000 we 
have, for practical reasons, less than one tenth that depth of data. We want to 
choose a fixed multiple of mixing time as batch length for all of our batches. To 
decide what range of mixing times will be reliable we used various portions of the 
data from those of our longest runs, and then applied the conclusions we drew to 
the other 3/4 of the runs. 

More specifically, we treat as "reliable" the empirical volume fraction of the 
longest runs, those of length at least 800 times mixing time. We then consider a 
range of batch partitions of these systems to see which ones give accurate t-test 
results. We are looking for 95% confidence intervals, so we expect such intervals to 
contain the true volume fraction 95% of the time; since the true volume fraction 
is unknown we instead check how frequently the intervals contain the empirical 
volume fraction, which for the longer runs we have assumed is reliable. We do this 
for each of the runs of length 800 or more times mixing time. The results on these 
systems are the following. 

For each of our longer runs (of at least 800 mixing times), we considered various 
initial portions of the run in each of six ranges of mixing times: 20-100, 100-200, 
300-400, 400-500, and 500-600. For each of these truncated runs we considered 
batch partitions of the data into equal size batches of a variety of multiples of 
mixing time: 1-5, 6-10, 11-15, 16-20, 21-30, 31-40 and 41-50. For each size run and 
for each batch size we computed a 95% confidence interval for the true mean of the 
volume fraction, and determined whether or not the confidence interval covers the 
sample mean for the full run (which we are assuming is interchangeable with the 
true mean) . The fraction of the more than 200 cases in each category for which the 
sample mean lies within the confidence interval is recorded in Table 5. From this it 
appears that using batches of size 1-5 mixing times would be unreliable, but that 
size 10 times mixing times would be reliable. (Table 5 is based on mixing times 
obtained using the autocorrelation /2. Table 6 is similar, using the autocorrelation 
/i, and again justifies the use of batches of size 10 times mixing time.) 
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We then used batches of size 10 times mixing times to obtain 95% confidence 
intervals for the true mean of all the systems, obtaining the results tabulated in 
Table 3 and included in Figs. 17 and 18. 

Finally, we applied the above batch criterion to obtain 90% confidence intervals 
for the true standard deviation of all our systems, using the method described earlier 
in this section. The results are in Table 4 and in Figs. 7 to 10. 

4. Conclusion 

We have performed Markov chain Monte Carlo simulations on a two dimen- 
sional model of low pressure granular matter of the general Edwards probabilistic 
type [EO]. Our main result, superficially summarized in Fig. 8, is that in this model 
the standard deviation of the volume fraction decays to zero as the particle num- 
ber increases, which indicates a well-defined random loose packing density for the 
model. This suggests that real granular matter exhibits sharply defined random 
loose packing; this could be verified by repeating the sedimentation experiments 
[JS] at a range of physical dimensions. Our argument is only convincing to the 
extent that the confidence intervals in Fig. 8 are small and justified, which required 
a statistical treatment of the data unusual in the physics literature. We hope that 
our detailed error analysis may be useful in other contexts. 
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Figure 1. Plot of volume fraction versus number of moves, from an initial volume 
fraction of 0.9991, for 970 squares. 
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Figure 2. Plot of volume fraction versus number of moves, from an initial volume 
fraction of 0.5192, for 994 squares. 
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Figure 3. An allowed configuration 



Figure 4. A uniform configuration 
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Figure 5. 399 squares after 10^ moves. 
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Figure 6a. Plot of biased autocorrelation (for 8000 particles) versus number of 
moves, with horizontal lines denoting one standard deviation away from zero. 
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Figure 6b. Close-up of the above plot, with initialization time. 



14 




Figure 7. Plot of the standard deviation of the volume fraction versus number of 
squares, using /2 for confidence intervals. 
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Figure 8. Plot of the standard deviation of the volume fraction versus number of 
squares, using log scales and /2 for confidence intervals. The line is y = —0.5004 x — 
0.8052. 
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Figure 9. Plot of the standard deviation of the volume fraction versus number of 
squares, using /i for confidence intervals. 
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Figure 10. Plot of the standard deviation of the volume fraction versus number of 
squares, using log scales and /i for confidence intervals. The line is y = — 0.5003 x — 
0.8055. 
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Figure 11. Plot of the mean of the volume fraction versus number of squares, using 
/2 for confidence intervals. The curve is y = 0.7541 + 0.0998 x~^/^. 
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Figure 12. Plot of the mean of the volume fraction versus number of squares, using 
/i for confidence intervals. The curve is y = 0.7541 + 0.0998 x~^/^. 
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Table 1 



Basics 

(using unbiased autocorrelation /i) 



number of number of 
squares moves in 
in packing units of 10^ 



99 


0.4 


100 


0.4 


195 


0.4 


205 


0.4 


294 


0.4 


294 


0.4 


399 


2 


413 


2 


497 


2 


504 


2 


600 


2 


603 


2 


690 


2 


690 


2 


790 


2 


803 


2 


913 


2 


913 


2 


996 


12 


1001 


12 


1955 


12 


2980 


12 


3003 


12 


3933 


12 


4008 


12 


4995 


12 


5908 


12 


6030 


12 


7037 


12 


7161 


12 


8015 


12 


8991 


12 


8995 


12 



step size ki 
in units (init. time) 
of 10^ in units of 
moves step size 
using /i 



0.2 


1 


0.2 


1 


0.2 


1 


0.2 


1 


0.2 


5 


0.2 


4 


0.2 


10 


0.2 


9 


0.2 


17 


0.2 


17 


0.2 


12 


0.2 


23 


0.2 


14 


0.2 


20 


0.2 


43 


0.2 


14 


0.2 


34 


0.2 


75 


1 


13 


1 


16 


1 


38 


1 


44 


1 


51 


1 


64 


1 


99 


1 


193 


1 


163 


1 


143 


1 


223 


1 


222 


1 


132 


1 


283 


1 


632 



kM run length 

(mixing time) in units of 

in units of mixing time 

step size using /i 
using /i 



1 


1998 


1 


1998 


1 


1998 


1 


1998 


5 


398 


4 


498 


10 


998 


9 


1110 


17 


587 


17 


587 


12 


832 


23 


433 


14 


713 


20 


498 


43 


231 


14 


713 


34 


293 


73 


135 


13 


953 


16 


755 


38 


318 


44 


277 


53 


234 


63 


193 


75 


158 


174 


68 


96 


125 


143 


84 


261 


46 


181 


48 


120 


100 


287 


41 


631 


18 
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Table 2 



Basics 

(using biased autocorrelation /2) 



number of number of 
squares moves in 
in packing units of 10^ 



99 


0.4 


100 


0.4 


195 


0.4 


205 


0.4 


294 


0.4 


294 


0.4 


399 


2 


413 


2 


497 


2 


504 


2 


600 


2 


603 


2 


690 


2 


690 


2 


790 


2 


803 


2 


913 


2 


913 


2 


996 


12 


1001 


12 


1955 


12 


2980 


12 


3003 


12 


3933 


12 


4008 


12 


4995 


12 


5908 


12 


6030 


12 


7037 


12 


7161 


12 


8015 


12 


8991 


12 


8995 


12 



step size kj 
in units (init. time) 
of 10^ in units of 
moves step size 
using /2 



0.2 


1 


0.2 


3 


0.2 


1 


0.2 


1 


0.2 


5 


0.2 


4 


0.2 


12 


0.2 


11 


0.2 


20 


0.2 


19 


0.2 


12 


0.2 


23 


0.2 


14 


0.2 


25 


0.2 


43 


0.2 


15 


0.2 


36 


0.2 


77 


1 


18 


1 


16 


1 


38 


1 


56 


1 


51 


1 


76 


1 


101 


1 


191 


1 


179 


1 


139 


1 


220 


1 


199 


1 


131 


1 


282 


1 


565 



kM run length 

(mixing time) in units of 

in units of mixing time 

step size using /i 
using /2 



1 


1998 


2 


998 


1 


1998 


1 


1998 


5 


398 


4 


498 


12 


832 


11 


908 


20 


498 


19 


525 


12 


832 


23 


433 


15 


665 


20 


498 


44 


226 


14 


713 


38 


262 


75 


132 


18 


687 


16 


755 


38 


318 


56 


218 


54 


230 


78 


156 


88 


135 


177 


67 


100 


120 


151 


80 


290 


41 


186 


47 


126 


95 


355 


33 


609 


19 
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Table 3 



Volume Fraction 



number of 

squares 
in packing 



sample value 



end points of 
95% confidence 
interval for 

true value, 
using /i 



end points of 
95% confidence 
interval for 
true value, 
using /2 



99 


0, 


,7637 


0, 


,7637±0, 


.0007 


0.7637±0, 


.0007 


100 


0. 


,7631 


0. 


,7631±0, 


.0007 


0.7631±0, 


.0007 


195 


0. 


,7617 


0. 


,7617±0, 


.0005 


0.7617±0, 


.0005 


205 


0. 


,7608 


0. 


.7608±0, 


.0005 


0.7608±0. 


.0005 


294 


0. 


,7605 


0. 


,7605±0, 


.0004 


0.7605±0, 


,0004 


294 


0. 


,7598 


0. 


,7598±0, 


.0004 


0.7598±0. 


.0004 


399 


0. 


,7596 


0. 


,7596±0, 


.0002 


0.7596±0, 


.0002 


413 


0. 


,7593 


0. 


,7593±0, 


.0002 


0.7593±0, 


.0002 


497 


0, 


,7590 


0, 


,7590±0, 


.0002 


0.7590±0, 


.0002 


504 


0. 


,7590 


0. 


,7590±0, 


.0003 


0.7590±0. 


.0002 


600 


0. 


,7583 


0. 


,7583±0, 


.0002 


0.7583±0, 


.0002 


603 


0. 


,7588 


0. 


,7588±0, 


.0002 


0.7588±0, 


,0002 


690 


0. 


,7581 


0. 


,7581±0, 


.0002 


0.7581±0, 


,0002 


690 


0. 


,7583 


0. 


,7583±0, 


.0003 


0.7583±0, 


,0003 


790 


0. 


,7578 


0. 


,7578±0, 


.0003 


0.7578±0. 


.0003 


803 


0. 


,7579 


0. 


,7579±0, 


.0002 


0.7579±0, 


,0002 


913 


0. 


,7575 


0. 


,7575±0, 


.0002 


0.7575±0. 


.0003 


913 


0. 


,7578 


0. 


,7578±0, 


.0003 


0.7578±0. 


.0003 


996 


0. 


,7575 


0. 


,7575±0, 


.0001 


0.7575±0, 


.0001 


1001 


0. 


,7574 


0. 


,7574±0, 


.0001 


0.7574±0, 


,0001 


1955 


0. 


,7565 


0. 


,7565±0, 


.0002 


0.7565±0. 


,0002 


2980 


0. 


,7558 


0. 


,7558±0, 


.0002 


0.7558±0, 


.0001 


3003 


0. 


,7559 


0. 


,7559±0, 


.0002 


0.7559±0, 


,0003 


3933 


0. 


,7554 


0. 


,7554±0, 


.0002 


0.7554±0, 


,0002 


4008 


0. 


,7558 


0. 


,7558±0, 


.0003 


0.7558±0. 


,0003 


4995 


0. 


,7555 


0. 


,7555±0, 


.0004 


0.7555±0, 


.0005 


5908 


0. 


,7549 


0. 


,7549±0, 


.0003 


0.7549±0, 


,0003 


6030 


0. 


,7551 


0. 


,7551±0, 


.0004 


0.7551±0, 


.0003 


7037 


0. 


,7545 


0. 


,7545±0, 


.0005 


0.7545±0. 


.0004 


7161 


0. 


,7550 


0. 


,7550±0, 


.0005 


0.7550±0. 


.0007 


8015 


0. 


,7551 


0. 


,7551±0, 


.0004 


0.7551±0, 


.0004 


8991 


0. 


,7551 


0. 


,7551±0, 


.0004 


0.7551±0. 


.0012 


8995 


0. 


,7550 




none 




none 
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Table 4 



Standard deviation of volume fraction 



number of 

squares 
in packing 



sample value 



end points of 
90% confidence 
interval for 

true value, 
using /i 



end points of 
90% confidence 
interval for 
true value, 
using /2 



99 





.0160 


0, 


.0160±0, 


.0005 


o.oieozbo, 


.0005 


100 


0, 


.0157 


0, 


.0157±0, 


.0006 


0.0157±0, 


.0006 


195 


0, 


.0110 


0, 


.OllOiO, 


.0004 


O.OllOiO, 


.0004 


205 


0, 


.0108 


0, 


.0108±0, 


.0004 


0.0108±0. 


.0004 


294 





.0090 


0, 


.0090±0, 


.0003 


0.0090±0, 


,0003 


294 


0, 


.0090 


0, 


.0090±0, 


.0003 


0.0090±0. 


.0003 


399 


0, 


.0078 


0, 


.0078±0, 


.0001 


0.0078±0, 


.0002 


413 





.0076 





.0077±0, 


.0001 


0.0077±0, 


.0001 


497 





.0070 





.0070±0, 


.0001 


0.0070±0, 


.0001 


504 


0, 


.0069 


0, 


.0069±0, 


.0001 


0.0069±0. 


.0001 


600 


0, 


.0064 


0, 


.0064±0, 


.0001 


0.0064±0, 


.0001 


603 


0, 


.0064 





.0064±0, 


.0001 


0.0064±0, 


,0001 


690 


0, 


.0060 


0, 


.0061±0, 


.0002 


0.0060±0, 


,0001 


690 


0, 


.0060 


0, 


.0060±0, 


.0002 


o.ooeozbo, 


,0001 


790 


0, 


.0055 


0, 


.0055±0, 


.0002 


0.0055±0, 


.0002 


803 





.0055 





.0055±0, 


.0002 


0.0055±0, 


,0002 


913 


0, 


.0052 


0, 


.0052±0, 


.0002 


0.0052±0. 


.0002 


913 


0, 


.0053 


0, 


.0053±0, 


.0001 


0.0053±0. 


.0001 


996 


0, 


.0050 


0, 


.0050±0, 


.0001 


0.0050±0. 


.0001 


1001 





.0049 





.0049±0, 


.0001 


0.0049±0, 


,0001 


1955 


0, 


.0036 


0, 


.0036±0, 


.0001 


0.0036±0, 


.0001 


2980 


0, 


.0028 


0, 


.0028±0, 


.0001 


0.0028±0, 


.0001 


3003 


0, 


.0029 


0, 


.0029±0, 


.0001 


0.0029±0, 


,0001 


3933 


0, 


.0024 


0, 


.0025±0, 


.0001 


0.0025±0, 


,0001 


4008 


0, 


.0025 


0, 


.0025±0, 


.0001 


0.0025±0. 


.0002 


4995 


0, 


.0022 


0, 


.0022±0, 


.0002 


0.0022±0, 


.0002 


5908 


0, 


.0021 


0, 


.0021±0, 


.0002 


0.0021±0, 


,0002 


6030 


0, 


.0020 


0, 


.0020±0, 


.0002 


0.0020±0, 


.0002 


7037 


0, 


.0018 


0, 


.0019±0, 


.0002 


0.0019±0. 


.0002 


7161 


0, 


.0019 


0, 


.0020±0, 


.0003 


0.0021±0. 


.0004 


8015 


0, 


.0016 


0, 


.0017±0, 


.0002 


0.0017±0, 


,0002 


8991 


0, 


.0017 


0, 


.0017±0, 


.0003 


0.0021±0. 


.0008 


8995 


0, 


.0016 




none 




none 
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Table 5 

Fraction of times the given batch size gives acceptable confidence interval for 
given segment of total data of long runs, using unbiased autocorrelation /i 







20-100 

mixing 
times of 
total data 


100-200 

mixing 
times of 
total data 


200-300 

mixing 
times of 
total data 


300-400 

mixing 
times of 
total data 


400-500 

mixing 
times of 
total data 


500-600 

mixing 
times of 
total dat 




1-5 


0.0849 


0.0867 


0.1119 


0.1191 


0.1938 


0.2082 


number of 


6-10 


0.9410 


0.9394 


0.9830 


1.0000 


1.0000 


1.0000 


mixing 


11-15 


0.9648 


0.9231 


0.9656 


1.0000 


1.0000 


1.0000 


times 


16-20 


0.9524 


0.9095 


0.9777 


0.9879 


1.0000 


1.0000 


per batch 


21-31 


0.9650 


0.9177 


0.9673 


1.0000 


1.0000 


1.0000 




31-40 


0.9712 


0.9042 


0.9643 


1.0000 


1.0000 


0.9957 




41-51 


0.9375 


0.8869 


0.9402 


0.9511 


0.9783 


0.9402 



Table 6 

Fraction of times the given batch size gives acceptable confidence interval for 
given segment of total data of long runs, using biased autocorrelation /2 







20-100 
mixing 

times of 
total data 


100-200 
mixing 

times of 
total data 


200-300 
mixing 

times of 
total data 


300-400 
mixing 

times of 
total data 


400-500 
mixing 

times of 
total data 


500-600 
mixing 

times of 
total datf 




1-5 


0.0833 


0.0924 


0.1182 


0.1259 


0.2006 


0.2326 


number of 


6-10 


0.9245 


0.9784 


0.9805 


0.9552 


0.9762 


1.0000 


mixing 


11-15 


0.9107 


0.9451 


0.9607 


0.9524 


0.9707 


1.0000 


times 


16-20 


0.9728 


0.9212 


0.9745 


0.9493 


0.9655 


1.0000 


per batch 


21-31 


0.9486 


0.9059 


0.9592 


0.9547 


0.9744 


1.0000 




31-40 


0.9780 


0.9190 


0.9592 


0.9704 


0.9852 


1.0000 




41-51 


0.9000 


0.8639 


0.9441 


0.9565 


0.9814 


1.0000 
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